title: “Informe PISA” author: “Laura Martínez González de Aledo” date: “11/7/2020” output: html_document
El siguiente estudio correspone a un analisis del Informe PISA del año 2006. El conjunto de datos se ha construido utilizando la puntuación media en Ciencias por país del Programa para la Evaluación Internacional de Estudiantes, junto con el GNI per cápita (paridad del poder adquisitivo, dólares de 2005), el índice educativo, el índice de salud y el índice de desarrollo humano de la ONU (HDI).
El objetivo es modelizar la relación entre la puntuación media (OSS) y el resto de variables, utilizando modelos de splines y GAM. Se debe realizar CV cuando se pueda.
library(tidyverse)
library(broom) # modelos en df
library(flextable) # Tablas formateadas
library(mgcv) # para estimar gam
library(reshape2) # melt
library(magrittr) # Pipe operators
library(janitor) # Clean names
library(skimr) # == Summarize
library(imputeTS) # para cambiar los valores nulos por la media
library(PerformanceAnalytics) # grafico de correlaciones
library(corrplot) # para ver outliers
library(gam) # para calcular el gam
library(rsample) # para el split
library(boot) # cv para modelos glm()
En primer lugar importo la tabla, a continuación elimino las columnas no voy a utilizar. Y por último hago una limpieza de tabla, borrando los duplicados y cambiando los Na´s por la media en vez de eliminarlos porque perderíamos información.
pisa2006 <- read.csv("pisasci2006.csv")
pisa2006 %<>% clean_names()
colnames(pisa2006) # Ahora los nombres de las columnas esta en minuscula
## [1] "country" "overall" "issues" "explain" "evidence" "interest"
## [7] "support" "income" "health" "edu" "hdi"
pisa2006 %<>% distinct(country,.keep_all= TRUE) # duplicados
summarise_all(pisa2006, funs(sum(is.na(.)))) # contamos valores nulos
## country overall issues explain evidence interest support income health edu
## 1 0 8 8 8 8 8 8 4 4 6
## hdi
## 1 6
pisa2006 <- na_mean(pisa2006) # cambiamos los NaN por la media
datos_pisa <- select(pisa2006, -issues, -explain, -evidence)
# quitamos las variables categoricas
head(datos_pisa) # vemos que nuestra no varia porque no teniamos duplicados
## country overall interest support income health edu hdi
## 1 Albania 473.1404 528.1579 512.1754 0.599 0.886 0.716000 0.7240000
## 2 Argentina 391.0000 567.0000 506.0000 0.678 0.868 0.786000 0.7730000
## 3 Australia 527.0000 465.0000 487.0000 0.826 0.965 0.978000 0.9200000
## 4 Austria 511.0000 507.0000 515.0000 0.835 0.944 0.824000 0.8660000
## 5 Azerbaijan 382.0000 612.0000 542.0000 0.566 0.780 0.799322 0.8060169
## 6 Belgium 510.0000 503.0000 492.0000 0.831 0.935 0.868000 0.8770000
# libreria skimr
skim(datos_pisa)
| Name | datos_pisa |
| Number of rows | 65 |
| Number of columns | 8 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| numeric | 7 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| country | 0 | 1 | 4 | 24 | 0 | 65 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| overall | 0 | 1 | 473.14 | 51.06 | 322.00 | 436.00 | 487.00 | 511.00 | 563.00 | ▁▂▂▇▃ |
| interest | 0 | 1 | 528.16 | 46.62 | 448.00 | 503.00 | 528.16 | 549.00 | 644.00 | ▆▇▇▃▂ |
| support | 0 | 1 | 512.18 | 24.40 | 447.00 | 497.00 | 512.18 | 527.00 | 569.00 | ▂▃▇▅▂ |
| income | 0 | 1 | 0.74 | 0.11 | 0.41 | 0.66 | 0.74 | 0.83 | 0.94 | ▁▂▆▇▅ |
| health | 0 | 1 | 0.89 | 0.06 | 0.72 | 0.84 | 0.89 | 0.94 | 0.99 | ▂▂▇▇▇ |
| edu | 0 | 1 | 0.80 | 0.11 | 0.54 | 0.74 | 0.80 | 0.87 | 0.99 | ▂▃▇▇▅ |
| hdi | 0 | 1 | 0.81 | 0.08 | 0.58 | 0.75 | 0.81 | 0.87 | 0.94 | ▁▃▆▇▆ |
Observamos las correlaciones para ver si hay problemas de multicolinealidad, dependiendo de si su correlación es lineal o no. Quito la variable Country porque no es una variable númerica
# librería corrplot
corrplot(cor(datos_pisa%>%
select_at(vars(-country)),
use = "complete.obs"),
method = "circle", type = "upper")
# libreria PerformanceAnalytics
chart.Correlation(datos_pisa %>%
select_at(vars(-country)),
histogram=TRUE, pch=12)
# Grafico interest
baseplot1 <- ggplot(data = datos_pisa, mapping = aes(x = overall, y = interest)) +
layer(geom = "point",stat = "identity",position = "identity") +
theme_bw() + theme(legend.key = element_blank())
baseplot1
# Grafico Support
baseplot1 <- ggplot(data = datos_pisa, mapping = aes(x = overall, y = support)) +
layer(geom = "point",stat = "identity",position = "identity") +
theme_bw() + theme(legend.key = element_blank())
baseplot1
# Grafico income
baseplot1 <- ggplot(data = datos_pisa, mapping = aes(x = overall, y = income)) +
layer(geom = "point",stat = "identity",position = "identity") +
theme_bw() + theme(legend.key = element_blank())
baseplot1
# Grafico health
baseplot1 <- ggplot(data = datos_pisa, mapping = aes(x = overall, y = health)) +
layer(geom = "point",stat = "identity",position = "identity") +
theme_bw() + theme(legend.key = element_blank())
baseplot1
# Grafico edu
baseplot1 <- ggplot(data = datos_pisa, mapping = aes(x = overall, y = edu)) +
layer(geom = "point",stat = "identity",position = "identity") +
theme_bw() + theme(legend.key = element_blank())
baseplot1
# Grafico hdi
baseplot1 <- ggplot(data = datos_pisa, mapping = aes(x = overall, y = hdi)) +
layer(geom = "point",stat = "identity",position = "identity") +
theme_bw() + theme(legend.key = element_blank())
baseplot1
A continuación calculamos los grados de libertad de cada variable junto con el Cross Validation. Unicamente los calculo para las variables númericas.
datos_pisa <- select(pisa2006, -country, -issues, -explain, -evidence)
La función smooth.spline() de R permite ajustar smooth splines de forma sencilla, con la ventaja añadida de que el valor óptimo de smothness (λ) puede identificarse por cross-validation.
plot(datos_pisa$interest, datos_pisa$overall, col='gray')
fit1 <- smooth.spline(datos_pisa$interest, datos_pisa$overall, df=16)
fit2 <- smooth.spline(datos_pisa$interest, datos_pisa$overall, cv=TRUE)
lines(fit1, col='red', lwd=2)
lines(fit2, col='blue', lwd=1)
fit2$df # 4.75
## [1] 4.750171
plot(datos_pisa$support, datos_pisa$overall, col='gray')
fit3 <- smooth.spline(datos_pisa$support, datos_pisa$overall, df=16)
fit4 <- smooth.spline(datos_pisa$support, datos_pisa$overall, cv=TRUE)
lines(fit3, col='red', lwd=2)
lines(fit4, col='blue', lwd=1)
fit4$df # 2.001
## [1] 2.001243
plot(datos_pisa$income, datos_pisa$overall, col='gray')
fit5 <- smooth.spline(datos_pisa$income, datos_pisa$overall, df=16)
fit6 <- smooth.spline(datos_pisa$income, datos_pisa$overall, cv=TRUE)
lines(fit5, col='red', lwd=2)
lines(fit6, col='blue', lwd=1)
fit6$df # 4.24
## [1] 4.244952
plot(datos_pisa$health, datos_pisa$overall, col='gray')
fit7 <- smooth.spline(datos_pisa$health, datos_pisa$overall, df=16)
fit8 <- smooth.spline(datos_pisa$health, datos_pisa$overall, cv=TRUE)
lines(fit7, col='red', lwd=2)
lines(fit8, col='blue', lwd=1)
fit8$df # 2.002
## [1] 2.002844
plot(datos_pisa$edu, datos_pisa$overall, col='gray')
fit9 <- smooth.spline(datos_pisa$edu, datos_pisa$overall, df=16)
fit10 <- smooth.spline(datos_pisa$edu, datos_pisa$overall, cv=TRUE)
lines(fit9, col='red', lwd=2)
lines(fit10, col='blue', lwd=1)
fit10$df # 2.00
## [1] 2.002385
plot(datos_pisa$hdi, datos_pisa$overall, col='gray')
fit11 <- smooth.spline(datos_pisa$hdi, datos_pisa$overall, df=16)
fit12 <- smooth.spline(datos_pisa$hdi, datos_pisa$overall, cv=TRUE)
lines(fit11, col='red', lwd=2)
lines(fit12, col='blue', lwd=1)
fit12$df # 8.603
## [1] 8.603228
# suavizado de splines
fit2$lambda
## [1] 0.00465548
fit4$lambda
## [1] 81.82485
fit6$lambda
## [1] 0.007177517
fit8$lambda
## [1] 44.19162
fit10$lambda
## [1] 52.19016
fit12$lambda
## [1] 0.0002060338
# library(gam)
gam1 <- gam(overall~ s(interest, 4.75) + s(support, 2.001) + s(income, 4.24) + s(health, 2.002) + s(edu, 2.002) + s(hdi, 8.603), data=datos_pisa)
summary(gam1)
##
## Call: gam(formula = overall ~ s(interest, 4.75) + s(support, 2.001) +
## s(income, 4.24) + s(health, 2.002) + s(edu, 2.002) + s(hdi,
## 8.603), data = datos_pisa)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -60.486 -8.202 3.138 11.680 51.442
##
## (Dispersion Parameter for gaussian family taken to be 596.3)
##
## Null Deviance: 166828.9 on 64 degrees of freedom
## Residual Deviance: 24091.7 on 40.402 degrees of freedom
## AIC: 620.1483
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## s(interest, 4.75) 1.000 90398 90398 151.5983 2.956e-15 ***
## s(support, 2.001) 1.000 2921 2921 4.8977 0.03260 *
## s(income, 4.24) 1.000 4151 4151 6.9614 0.01177 *
## s(health, 2.002) 1.000 745 745 1.2494 0.27027
## s(edu, 2.002) 1.000 342 342 0.5731 0.45340
## s(hdi, 8.603) 1.000 3320 3320 5.5670 0.02322 *
## Residuals 40.402 24092 596
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Anova for Nonparametric Effects
## Npar Df Npar F Pr(F)
## (Intercept)
## s(interest, 4.75) 3.8 3.3660 0.020073 *
## s(support, 2.001) 1.0 1.6190 0.210521
## s(income, 4.24) 3.2 6.2341 0.001094 **
## s(health, 2.002) 1.0 1.4042 0.243012
## s(edu, 2.002) 1.0 1.7007 0.199599
## s(hdi, 8.603) 7.6 1.2307 0.307436
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(gam1, se=TRUE, col='blue')
gam2 <- gam(overall~ s(interest, 4.75) + s(support, 2.001)+ s(health, 2.002) + s(hdi, 8.603), data=datos_pisa)
summary(gam2)
##
## Call: gam(formula = overall ~ s(interest, 4.75) + s(support, 2.001) +
## s(health, 2.002) + s(hdi, 8.603), data = datos_pisa)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -81.927 -11.914 3.562 15.920 48.927
##
## (Dispersion Parameter for gaussian family taken to be 720.8222)
##
## Null Deviance: 166828.9 on 64 degrees of freedom
## Residual Deviance: 33622.08 on 46.6441 degrees of freedom
## AIC: 629.3297
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## s(interest, 4.75) 1.000 103182 103182 143.1452 8.117e-16 ***
## s(support, 2.001) 1.000 2016 2016 2.7975 0.10111
## s(health, 2.002) 1.000 4341 4341 6.0220 0.01792 *
## s(hdi, 8.603) 1.000 417 417 0.5788 0.45061
## Residuals 46.644 33622 721
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Anova for Nonparametric Effects
## Npar Df Npar F Pr(F)
## (Intercept)
## s(interest, 4.75) 3.8 3.7802 0.01094 *
## s(support, 2.001) 1.0 2.1137 0.15267
## s(health, 2.002) 1.0 0.9677 0.33050
## s(hdi, 8.603) 7.6 1.5333 0.17476
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(gam2, se=TRUE, col='blue')
gam3 <- gam(overall~ s(interest, 4.75) + s(support, 2.001) + s(edu, 2.002), data=datos_pisa)
summary(gam3)
##
## Call: gam(formula = overall ~ s(interest, 4.75) + s(support, 2.001) +
## s(edu, 2.002), data = datos_pisa)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -66.389 -10.251 5.461 12.100 64.652
##
## (Dispersion Parameter for gaussian family taken to be 743.8547)
##
## Null Deviance: 166828.9 on 64 degrees of freedom
## Residual Deviance: 41095.59 on 55.2468 degrees of freedom
## AIC: 625.1709
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## s(interest, 4.75) 1.000 103773 103773 139.5065 < 2e-16 ***
## s(support, 2.001) 1.000 2632 2632 3.5383 0.06524 .
## s(edu, 2.002) 1.000 2821 2821 3.7930 0.05656 .
## Residuals 55.247 41096 744
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Anova for Nonparametric Effects
## Npar Df Npar F Pr(F)
## (Intercept)
## s(interest, 4.75) 3.8 4.4461 0.004134 **
## s(support, 2.001) 1.0 2.6246 0.110885
## s(edu, 2.002) 1.0 1.4429 0.234867
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(gam3, se=TRUE, col='blue')
Realizamos el test anova para comparar los tres modelos que hemos propuesto anteriormente.
anova.gam(gam1, gam2, gam3, test='F')
## Analysis of Deviance Table
##
## Model 1: overall ~ s(interest, 4.75) + s(support, 2.001) + s(income, 4.24) +
## s(health, 2.002) + s(edu, 2.002) + s(hdi, 8.603)
## Model 2: overall ~ s(interest, 4.75) + s(support, 2.001) + s(health, 2.002) +
## s(hdi, 8.603)
## Model 3: overall ~ s(interest, 4.75) + s(support, 2.001) + s(edu, 2.002)
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 65 24092
## 2 65 33622 0 -9530.4
## 3 65 41096 0 -7473.5
Podemos comprobar que el que menor residuo tiene es el modelo 1 por lo que va a ser el modelo con el que vamos a trabajar.
Una vez escogido el modelo, vamos a dividir nuestra base de datos entre muestra para test y muestra para entrenamiento, vamos a proceder a introducir nuestro modelo en el test para saber como predice.
Dividimos aleatoriamente las observaciones disponibles en dos grupos, uno se emplea para entrenar al modelo y otro para evaluarlo.
# library(rsample)
set.seed(123)
pisa_split <- initial_split(datos_pisa, prop =.7, strata = "overall")
pisa_train <- training(pisa_split)
pisa_test <- testing(pisa_split)
Tenemos la base de datos dividida en 70/30, y vamos a proceder a introducir nuestro modelo en el test para saber como predice.
gam_train <- gam(overall~ s(interest, 4.75) + s(support, 2.001) + s(income, 4.24) + s(health, 2.002) + s(edu, 2.002) + s(hdi, 8.603), data=pisa_train)
plot(gam_train, se=TRUE, col='red')
summary(gam_train)
##
## Call: gam(formula = overall ~ s(interest, 4.75) + s(support, 2.001) +
## s(income, 4.24) + s(health, 2.002) + s(edu, 2.002) + s(hdi,
## 8.603), data = pisa_train)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -28.8039 -8.0293 0.5306 8.2329 27.3118
##
## (Dispersion Parameter for gaussian family taken to be 332.6423)
##
## Null Deviance: 108148.2 on 46 degrees of freedom
## Residual Deviance: 7451.927 on 22.4022 degrees of freedom
## AIC: 422.6816
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## s(interest, 4.75) 1.000 59311 59311 178.3015 3.793e-12 ***
## s(support, 2.001) 1.000 1418 1418 4.2632 0.05071 .
## s(income, 4.24) 1.000 8218 8218 24.7065 5.387e-05 ***
## s(health, 2.002) 1.000 1921 1921 5.7749 0.02496 *
## s(edu, 2.002) 1.000 47 47 0.1422 0.70966
## s(hdi, 8.603) 1.000 287 287 0.8622 0.36302
## Residuals 22.402 7452 333
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Anova for Nonparametric Effects
## Npar Df Npar F Pr(F)
## (Intercept)
## s(interest, 4.75) 3.8 3.1114 0.037830 *
## s(support, 2.001) 1.0 0.9666 0.336116
## s(income, 4.24) 3.2 6.1667 0.002732 **
## s(health, 2.002) 1.0 5.2212 0.032069 *
## s(edu, 2.002) 1.0 0.1249 0.727573
## s(hdi, 8.603) 7.6 2.5807 0.038446 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
La función cv.glm() calcula el error de predicción mediante cross-validation.
# library(boot)
cv <- cv.glm(data=datos_pisa, glmfit = gam1)
sqrt(cv$delta)
## [1] 35.19908 35.05841
Igual que el modelo anterior sólo que especificando K.
set.seed(123)
cv_error <- cv.glm(data=datos_pisa, glmfit = gam1, K=10)
cv_error$delta
## [1] 1250.154 1189.006
Cuando se especifica el número de grupos K a emplear en la validación, la función cv.glm() devuelve dos resultados, uno con corrección de continuidad y otro sin ella.
Smoothing Splines: https://rpubs.com/Joaquin_AR/250069